library(sasld)


# 10.1 ---------------------------------------------------------------
# Esempio 4 (Tabella 10.2)
prop.test(c(347,327),c(11535,14035),correct=FALSE)

tmp <- prop.test(c(347,327),c(11535,14035),correct=FALSE)
tmp$conf.int
round(tmp$conf.int,3)

tmp <- prop.test(c(347,327),c(11535,14035),conf.level=0.99,correct=FALSE)
round(tmp$conf.int,3)

num <- c(347,327)
den <- c(11535,14035)
p <- num/den
p[1] - p[2]
varp <- p*(1-p)/den
sqrt(varp[1] + varp[2])
# -------------------------------------------------------------------




# 10.2 ---------------------------------------------------------------
# Esempio 5 (Tabella 10.3)
prop.test(c(5,154),c(88,619),correct=FALSE)

# SE
num <- c(5,154)
den <- c(88,619)
p <- num/den
(ph <- sum(num)/sum(den))
vph <- ph*(1-ph)/den
(se <- sqrt(sum(vph)))
# -------------------------------------------------------------------


# 10.3 ---------------------------------------------------------------
# simulazione (pag. 473 testo inglese)
set.seed(123456)
nrep <- 100000
x1 <- rbinom(n=nrep,size=11535,prob=0.0265)
x2 <- rbinom(n=nrep,size=14035,prob=0.0265)
p1 <- x1/11535
p2 <- x2/14035

cbind(mean(p1),mean(p2),mean((p1-p2)))
cbind(sd(p1),sd(p2))
#
sqrt(0.0265*(1-0.0265)/11535)
sqrt(0.0265*(1-0.0265)/14035)

sd((p1-p2))
#
sqrt((0.0265*(1-0.0265)/11535) + (0.0265*(1-0.0265)/14035))

err1 <- p1 - 0.0265
err2 <- p2 - 0.0265
delta <- p1 - p2

cbind(mean(err1),mean(err2),mean(delta))

diff.1 <- abs(delta) - abs(err1)
diff.2 <- abs(delta) - abs(err2)

mean(diff.1)
mean(diff.2)

table(sign(diff.1))
table(sign(diff.2))
# -------------------------------------------------------------------


# 10.4 ---------------------------------------------------------------
tmp <- prop.test(c(5+1,154+1),c(88+2,619+2),correct=FALSE)
round(tmp$conf.int,3)

tmp <- prop.test(c(1,1),c(12,12),correct=FALSE)
round(tmp$conf.int,3)
# -------------------------------------------------------------------


# 10.5 ---------------------------------------------------------------
# Esempio 8
test.t2ci(media1=5.9,ds1=3.3,n1=75,media2=1,ds2=2.3,n2=257)
# -------------------------------------------------------------------


# 10.6 ---------------------------------------------------------------
# Esempio 9
test.t2ci(media1=585.2,ds1=89.6,n1=32,media2=533.7,ds2=65.3,n2=32)

# escludendo l'outlier
test.t2ci(media1=573.1,ds1=58.9,n1=31,media2=533.7,ds2=65.3,n2=32)
pt(2.516338,60.69685,lower.tail=FALSE) # nel testo è 0.008
# -------------------------------------------------------------------


# 10.7 ---------------------------------------------------------------
# Esempio 10
test.t2ci(media1=51.6,ds1=23.7,n1=60,media2=53.7,ds2=23.7,n2=61)
# -------------------------------------------------------------------


# 10.8 ---------------------------------------------------------------
# Esempio 13
data(drt)

t.test(drt$Yes,drt$No,paired=TRUE)
test.t2cd(cbind(drt$Yes,drt$No))

delta <- drt$Yes - drt$No
m <- mean(delta)
d <- sd(delta)
n <- length(delta)
test <- m/(d/sqrt(n))
c(m,d,n,test)
# -------------------------------------------------------------------


# 10.9 ---------------------------------------------------------------
# Esempio 15
a1 <- rep(1,955); b1 <- a1
a2 <- rep(1,162); b2 <- rep(0,162)
a3 <- rep(0,9);   b3 <- rep(1,9)
a4 <- rep(0,188); b4 <- a4

heaven <- c(a1,a2,a3,a4)
hell   <- c(b1,b2,b3,b4)
table(heaven,hell)

test.t2cd(cbind(heaven,hell))

# test di McNemar
tbl <- c(188,9,162,955)
tbl <- matrix(tbl,nrow=2,byrow=TRUE)
mcnemar.test(tbl, correct=FALSE)

((162-9)/sqrt(162+9))^2
# -------------------------------------------------------------------


# 10.10 -------------------------------------------------------------
# Esempio 16 (test esatto)
binom.test(16,74)

tbl <- c(1921,58,16,5)
tbl <- matrix(tbl,nrow=2,byrow=TRUE)
mcnemar.test(tbl, correct=FALSE)
# -------------------------------------------------------------------